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Abstract 

In order to perform numerical studies of long-term stability in nonlinear 
Hamiltonian systems, one needs a numerical integration algorithm which is 
symplectic. Further, this algorithm should be fast and accurate. In this paper, 
we propose such a symplectic integration algorithm using polynomial map 
refactorization of the symplectic map representing the Hamiltonian system. 
This method should be particularly useful in long-term stability studies of 
particle storage rings in accelerators. 
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I. INTRODUCTION 



Long-term single particle stability studies of particle storage rings play an important role 
in the design of accelerators [Q. These storage rings are generally described by nonlinear, 
nonintegrable Hamiltonians. Therefore analytical results on long-term stability of particle 
motion in such storage rings are difficult to obtain. By default, numerical integration of par- 
ticle trajectories is the primary tool used to explore the dynamics of these systems. However, 
standard numerical integration algorithms can not be used since they are not symplectic [Q. 
This violation of the symplectic condition can lead to spurious chaotic or dissipative be- 
havior. Numerical integration algorithms which satisfy the symplectic condition are called 
symplectic integration algorithms 0. 

Several symplectic integration algorithms have been proposed in the literature [^]. Some 
of these directly use the Hamiltonian whereas others use the symplectic map p[ representing 
either the entire storage ring (in which case one obtains the so-called one-turn map) or major 
segments of the ring. For complicated systems like the Large Hadron Collider which has 
thousands of elements, using individual Hamiltonians for each element can drastically slow 
down the integration process. One the other hand, the map based approach is very fast in 
such cases . Further, if nonlinearities in the symplectic map are too "large" , one can use 
scaling and squaring techniques to alleviate the problem. 

One class of the map-based methods uses jolt factorization 0J^. But there are still 
unanswered questions on how to best choose the underlying group and elements in the 
group 0. Further, some of these methods can be quite difficult to generalize to higher 
dimensions. Another class of methods uses solvable maps 0] or monomial maps HIT 



Even though they are fairly straightforward to generalize to higher dimensions, they tend 
to introduce spurious poles and branch points not present in the original map 0. 

In this paper, we propose a new symplectic integration method where the symplectic map 
is refactorized using "polynomial maps" (maps whose action on phase space variables gives 
rise to polynomials). This does not introduce spurious poles and branch points. Moreover, 
it is easy to generalize to higher dimensions. In this letter, we restrict ourselves to maps in 
six dimensional phase space which are appropriate for single particle stability studies. We 
show that the method gives good results. Further, since it is map-based, it is also very fast. 



II. PRELIMINARIES 

We start by representing a Hamiltonian system by a symplectic map 0. For simplicity 
we restrict ourselves to a six dimensional phase space. Let us denote the collection of six 
phase-space variables qi, pi {i = 1, 2, 3) by the symbol z: 

^ = (9l,?2,g3,Pl,P2,P3)- (1) 

The Lie operator corresponding to a phase-space function f{z) is denoted by -.f^z):. It 
is defined by its action on a phase-space function g{z) as shown below 

■.f{z):g{z) = [f{z),g{z)]. (2) 
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Here [f{z),g{z)] denotes the usual Poisson bracket of the functions f{z) and g{z). Next, we 
define the exponential of a Lie operator. It is called a Lie transformation and is given as 
follows: 

^.n.y. ^ £ iMi!. (3) 

n=0 "'• 

The effect of a Hamiltonian system on a particle can be formally expressed as the action 
of a map Ai that takes the particle from its initial state 2;™ to its final state z-^^"' 

z^''' = Mz'''. (4) 

It can be shown that is a symplectic map 0. Consider its Jacobian matrix which we 
denote by M. Symplectic maps are maps whose Jacobian matrices M satisfy the following 
'symplectic condition' 

MJM = J, (5) 

where M is the transpose of M and J is the fundamental symplectic matrix. 

Using the Dragt-Finn factorization theorem |T^fl, the symplectic map Ai can be fac- 
torized as shown below: 

M = Me-^^'- e-^^'- . . . e^"'- .... (6) 

Here fn{z) denotes a homogeneous polynomial (in z) of degree n uniquely determined by 
the factorization theorem. Further M gives the linear part of the map and hence has an 
equivalent representation in terms of the Jacobian matrix M of the map Ai 0: 

Mzi = M,jz, = {Mz)i. (7) 

The infinite product of Lie transformations exp(:/„:) {n = 3,4, . . .) in Eq. represents 
the nonlinear part of Ai. 

Using the above procedure, one can represent each element in the storage ring by a 
symplectic map. By concatenating these maps together, we obtain the so-called 'one- 
turn' map representing the entire storage ring. This concatenation is made possible by the 
Campbell-Baker-Hausdorff (CBH) theorem |T^. The one-turn map gives the final state 2;^^^ 



of a particle after one turn around the ring as a function of its initial state z^^^ : 

zW = >f^{o). (8) 

To obtain the state of a particle after turns, one has to merely iterate the above mapping 
times i.e. 



z 



W=7W"2(o). (9) 



Since Ai is explicitly symplectic, this gives a symplectic integration algorithm. Further, 
since the entire ring can be represented by a single (or at most a few) symplectic map(s), 
numerical integration of particle trajectories using symplectic maps is very fast. 
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III. SYMPLECTIC INTEGRATION USING POLYNOMIAL MAPS 



It is obvious that one can not use Ad in the form given in Eq. (|^) for any practical 
computations. It involves an infinite number of Lie transformations. Therefore, we have to 
truncate Ai by stopping after a finite number of Lie transformations: 

M ^ Me-^'^'- e-f^'- . ..e^^'-. (10) 

However, each exponential e''^"' in Ai still contains an infinite number of terms in its Taylor 
series expansion. One possible solution is to truncate the Taylor series generated by the Lie 
transformations to order P. But this violates the symplectic condition. 

We get around the above problem by refactorizing A4 in terms of simpler symplectic 
maps which can be evaluated exactly without truncation. We use 'polynomial maps' which 
give rise to polynomials when acting on the phase space variables. This avoids the problem 
of spurious poles and branch points present in generating function methods ||^, solvable map 
jlOl and monomial map |ll| refactorizations. To determine which symplectic maps give rise 



to polynomial mappings, consider exp(: h{z) :) where h{z) is a polynomial. Since all Lie 
transformations are symplectic maps 0, this is a symplectic map. Its action on phase space 
variables is equivalent to solving the Hamilton's equations of motion from time t = to 
t = — 1 using h{z) as the Hamiltonian. For example, consider the action of exp(: qf :) on qi, 
Pi in a two dimensional phase space. We first solve the following Hamilton's equations of 
motion: 



dqi dh 
dt dpi ' 
dpi dh 



(11) 



dt dqi ' 

with h = qf. Solving these simple equations, we obtain: 

qi{t) = qi{0), Pi{t)=pm-3qi{0)H, (12) 

where gi(0) and pi(0) denote the values of qi and pi at time t = 0. To obtain the action 
of the map exp(: qf :) on the phase space variables, set t = — 1 in the above equations and 
denote gi(— 1), pi(— 1) by gf*", pi*" and gi(0), pi(— 0) by g^", p^" respectively. Thus we get 

gf" = gr, pt=pT+mn'- (13) 

Using Eq. (^, we can easily verify that the above result is indeed correct. 

Using the above procedure, we can identify symplectic maps exp(: h{z) :) which give 
rise to polynomial mappings of the phase space variables into themselves. These results 



14 1 can be codified as the following simple principles which are easily generalized to higher 



dimensions also. 

1. All polynomials of the form h{z) where both a phase space variable and its canonically 
conjugate variable [|l^] do not occur simultaneously give rise to polynomial symplectic 
maps via exp(: h{z) :). 
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2. If a canonically conjugate pair {qi,Pi} is present in the polynomial h{z), it only appears 
either in the form a{z')qi + g{pi, z') or a{z')pi + g{qi, z') or its integer powers. Here z' 
denotes the collection of phase space variables {qj,Pk} with j ^ k ^ i. Further, a and 
g are polynomials in the indicated variables. 

We now return to the problem of symplectic integration. For the present, we restrict 
ourselves to one-turn symplectic maps in a two dimensional phase space truncated at order 4. 
The results obtained below can be generalized to higher orders using symbolic manipulation 
programs. The Dragt-Finn factorization of the symplectic map is given by: 

M = Me-f'''- e-f^\ (14) 

where 

/s = a.iql + a2qlPi + amvl + 04^?, 

/4 = a^q{ + a^qlpi + a^qlpl + a^qipl + a<ip\. (15) 

Here the coefficients ai, . . . , og can be explicitly computed given a Hamiltonian system p|. 
The above map captures the leading order nonlinearities of the system. Since the action 
of the linear part M on phase space variables is well known [cf. Eq. (|^)] and is already a 
polynomial action, we only refactorize the nonlinear part of map using polynomial maps. It 
turns out that we require 7 polynomial maps for this purpose: 

M^V = Me-^^- e-^^- ■■■e^'\ (16) 

where the numeral appearing in the subscript indexes the polynomial maps. The /ij's are 
given as follows: 

hi = biqf + hr^ql, /ig = + hp\, 

hs = (&2 + h)iqi+pif, /i4 = ih - h2){qi-pi)\ (17) 
h-, = {qi + Pi + h^plf , h& = {-qi - pi + Kqlf , hj = br{qi + pi)^. 

Here 6j's are at present unknown coefficients. By forcing the refactorized form V to equal 
the original map Ai up to order 4 and using the CBH theorem we can easily compute 
these unknown coefficients in terms of the known a^'s. These expressions are given in the 
Appendix. The explicit actions of the polynomial maps on the phase space variables are 
also given there. This completely determines the refactorized map V. Each exp(: hi :) is 
a polynomial map which can be evaluated exactly and is explicitly symplectic. Thus by 
using V instead of Ai in Eq. {^), we obtain an explicitly symplectic integration algorithm. 
Further, it is fast to evaluate and does not introduce spurious poles and branch points. The 
above factorization is not unique. However, the principles outlined earlier impose restrictions 
on the possible forms and this eases considerably the task of refactorization. Moreover, we 
require the coefficients bi to be polynomials in the known coefficients Oj. Otherwise this can 
lead to divergences when a^'s take on certain special values. Finally, we minimize the number 
of polynomial maps in the refactorized form. Our studies show that different polynomial 
map refactorizations obeying the above restrictions do not lead to any significant differences 
in their behavior. 
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The above refactorization has also been extended to symplectic maps in a six dimensional 
phase space truncated at order 4. In this case, we require 23 polynomial maps in the 
refactorization to make V equal to M. up to order 4: 

M^V = Me-^'^- e-^^- ■■■e-^^''-, (18) 

where the numeral appearing in the subscript indexes the polynomial maps. Since listing 
out the explicit forms of the /ij's and their coefficients is not particularly illuminating, we do 
not list them here. However, a FORTRAN program which implements the above polynomial 
map refactorization is available from the author upon request. 

We now analyze the leading order error committed in our method. In our method, we 
first truncate the symplectic map to a given order and then refactorize it using a product of 
polynomial maps. Both these stages give rise to errors. When we truncate the symplectic 
map M. at the nth order, we obtain 

Mn = Mexp(: h ■■) exp(: U ■) ■ • •exp(: /„ :). (19) 

The leading term that has been omitted is exp(: :). From properties of Lie transforma- 
tions and Lie operators 0, we have 

exp(: fn+i ■.)z = z + [fn+i, z]^ , (20) 

where [, ] denotes the usual Poisson bracket. Now, gives terms of the form 0. 

Thus error due to truncation of the symplectic map is of order z^. 

Next, we refactorize the truncated symplectic map Ain as a product of k polynomial 
maps: 

Mn = Mexp(: hi :) exp(: /12 :) • • .exp(: hk :). (21) 

These polynomial maps are obtained by first using the CBH series to combine the Lie 
transformations and then comparing with the original symplectic map. Both these maps 
are made to agree up to order n. Therefore, the leading error term is again of the form 
exp(: /„_!_! :) giving rise to an error of order z^. 



IV. APPLICATIONS 

We now consider two applications of the above method. The first example is to find the 
region of stability of the following simple symplectic map: 

M = Mexp[:{qi+pif:], (22) 

where 

V — sm cos 6 J 

and 6 = ^. We chose this example since the exact action of the above map is known and 
hence the exact region of stability can also be determined. We found excellent agreement 
between results obtained using polynomial maps and the exact results. 
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We have also applied the method to more complicated Hamiltonian systems like particle 
storage rings. We studied an electron storage ring with radio frequency bunching cavities. 
The storage ring is composed of drifts, bending magnets, quadrupoles, sextupoles and RF 
cavities. The efficacy of our method is best revealed for such complicated Hamiltonian 
systems. Since there are many constituent elements (in storage rings like the Large Hadron 
Collider, there can be thousands of elements), numerical integration using Hamiltonians for 
each element is cumbersome and slow. On the other hand, a map based approach where 
one represents the entire storage ring in terms of a single map is much faster |Q . When this 
is combined with our polynomial map refactorization, one obtains a symplectic integration 
algorithm which is both fast and accurate and is ideally suited for such complex real life 
systems. The gs — ps phase plot for one million turns around the ring using our polynomial 
map method is given in Figure 1. In this case, qs and ps represent the deviations from 
the closed orbit time of flight and energy respectively. From theoretical considerations, 
we expect the so-called synchrotron oscillations in these variables. This manifests itself as 
ellipses in the phase space plot of and ps variables (just as the oscillations of the simple 
pendulum manifest themselves as ellipses in the coordinate-momentum phase space plot). 
In Figure 1, we observe the expected synchrotron oscillations. 



V. CONCLUSIONS 

To conclude, we have proposed a new symplectic integration algorithm based on poly- 
nomial map refactorization. This should be of help in studying long term stability of com- 
plicated accelerator systems and other Hamiltonian systems. 
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APPENDIX 



The coefficients 6j in Eq. ( p!7D can be easily determined using the CBH theorem []I3| . 
Their expressions in terms of the known coefficients of the symplectic map Ai [cf. Eq. 
(p!5|)] is given as follows: 



bi 


= ai — cis/S, 62 = CL2/G, 






h 


= a^/Q, 64 = 04 — 02/3, 






h 


= (2a6 - 2a7 + as + I86162 - 366^ - 


- 366163 


+ 366^ + 96164 + I86264 - 186364)/6, 


h 


= (— cte + — ctg ~ I86162 + 3662 


+ I86163 - 366^ - 96164 - I86264 + 186364)/4, (24) 




= (ae - 2ar + 2as + I86162 - 366^ - 


- I86163 


+ 366^ + 96164 + 366264 - 186364)/6, 




= as - 96162 - 96^ + 96^ - 860 - 67 






b9 


= ag - 96^ + 96^ + 96364 - 67 - 3bs 
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Note that the formulas have been sequenced in such a way that once a given bi is evaluated, 
it is used in the formulas for the fej's following it. 

The actions of the polynomial maps exp(: hi :) {i = 1,2, ...,7) on the phase space 
variables qi, pi are easily evaluated using the procedure outlined in the main text (see the 
discussion before Eq. (pJ]). We obtain the following results: 




gi, e-'^'-pi =pi + 3biql + Ab^qf, 

qi - Sb^pj - Abgpl, e^^'-pi = pi, 

qi - 3(62 + &3)(gi + Pi)\ e-'^'-'pi =pi + 3(62 + b^){q: 

qi + Kh-h2){qi-Pi)\ e-'^^'-pi = Pi + 3ibs - b2)iq_ 

qi - ci(l + 268^1 + &8C1), e^^'-pi = pi + ci, 

qi + C2, e-^^'-pi =pi- C2(l - 266^1 - beC2), 

qi - 467(51 + e-^^--p^ =Pi+ 467(51 + pif, 



+ 




(25) 



where 



ci = 3(gi + pi + b^plf, C2 = 3(-gi -pi + bepjf. 



(26) 
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FIGURE LEGENDS 

Figure 1: This figure sliows tlie — pa pliase space plot for one million turns around a 
storage ring using the polynomial map method (only every 1000th point is plotted). 
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